rm(list=ls(all=TRUE))

library(haven)
library(estimatr)
library(RColorBrewer)
library(ggplot2)
library(scales)
library(patchwork)
library(mediation)

# Experimental Results
d0 <- read_dta("victimization.dta", encoding = "latin1")
outcomes <- c(
  "nat","pat",
  "system",
  "perception","hawk")
y.labs <- c(
  "Victim Sentiments",
  "Patriotic Sentiments",
  "Support for Current System",
  "Perception of Intention", 
  "Support for Hawkish Policies"
)
y.lims <- list(c(-0.4, 0.6), c(-0.4, 0.6),c(-0.4, 0.6),
               c(-0.4, 0.6), c(-0.4,0.6), c(-0.4,0.6))
treat <- c("tr_victim","tr_liberal","tr_olympic")
tr.labs <- c("Tr: Victimizing", "Tr: Historical","Tr: Olympic")
covars <- c("ideology_nati","ideology_poli","ideology_econ",
            "age","age_sq","female","hschool","jcollege","college","factor(prov)", 
            "income", "substatus","knwl", "minority", "ccp")

mycol <- brewer.pal(7,"Set1")[c(1,2,3,4,5,7)]
mypch <- c(16, 17, 15, 18, 10, 12)
nvars <- length(y.labs)
adj <- c(-0.5,0,0.5)*0.2
for (var in outcomes) {
  d0[,var] <- scale(d0[,var])
}
d.high <- d0[which(d0$eduyr>=16),]
d.low <- d0[which(d0$eduyr<16),]
d.list <- list(d0, d.low, d.high)
sample.labs <- c("Full Sample","Below College","College or Above")

# Figure 3: Treatment Effects on Victim and Patriotic Sentiments
# estimates of victim sentiments
formula <- as.formula(paste0(outcomes[1],"~",
                             paste(treat, collapse = "+"),"+",paste(covars, collapse = "+")))
est <- array(NA, dim = c(length(treat), 5, 2))
out <- lm_robust(formula, data = d0, se_type = "HC1")
out2 <- lm_robust(formula, data = d0, se_type = "HC1", alpha = 0.1)
x0 <- c(1:length(treat))+1
est[,,1] <- cbind(out$coefficients[x0], out$conf.low[x0], out$conf.high[x0], out2$conf.low[x0], out2$conf.high[x0])

# estimates of patriotic sentiments
formula2 <- as.formula(paste0(outcomes[2],"~",
                              paste(treat, collapse = "+"),"+",paste(covars, collapse = "+")))
out3 <- lm_robust(formula2, data = d0, se_type = "HC1")
out4 <- lm_robust(formula2, data = d0, se_type = "HC1",alpha = 0.1)
x0 <- c(1:length(treat))+1
est[,,2] <- cbind(out3$coefficients[x0], out3$conf.low[x0], out3$conf.high[x0],out4$conf.low[x0], out4$conf.high[x0])

## plotting
pdf(paste0("graphs/","victim_pat.pdf"), width = 8, height = 4)
par(mar = c(3,3,3,1))
plot(1, type = "n", xlim = c(0.5, 2.5), 
     ylim = c(-0.4,0.6), axes = FALSE)
box(); axis(2); 
axis(1, at = 1:2, labels = c("Victim Sentiments", "Patriotic Sentiments"), line = -.2, tick = FALSE)
abline(h = seq(-1,1,0.1), col = "#AAAAAA50") # treatment line color
abline(v = seq(0,4,0.1), col = "#AAAAAA50")
abline(h = 0, lwd = 5, col = "#DDDDDD")
for (i in 1:2) {
  for (j in 1:length(treat)) {
    x0 <- i + adj[j]
    lines(c(x0,x0), c(est[j,2,i], est[j,3,i]), lwd = 1, col = "grey30")
    lines(c(x0,x0), c(est[j,4,i], est[j,5,i]), lwd = 3, col = "grey30")
    points(x0, est[j,1,i], col = "grey30", pch = mypch[j], cex = 1.5)
  }
}
legend("top",legend=tr.labs, border=NA,bty="n",ncol=3, col="grey30",
       lwd=2,pch=mypch,cex=1,seg.len=2)
graphics.off()


# Figure 4 Treatment Effects on Victim and Patriotic Sentiments by Education
d.high <- d0[which(d0$eduyr>=16),]
d.low <- d0[which(d0$eduyr<16),]
d.list <- list(d.low, d.high)
sample.labs <- c("Below College","College or Above")

# estimate of victim and patriotic sentiments by education
for (i in 1:2) {
  d0[,outcomes[i]] <- scale(d0[,outcomes[i]])
  formula.m <- c(as.formula(paste0(outcomes[1],"~",
                                   paste(treat, collapse = "+"),"+",paste(covars, collapse = "+"))), 
                 as.formula(paste0(outcomes[2],"~",
                                   paste(treat, collapse = "+"),"+",paste(covars, collapse = "+"))))
  s <- d.list[[i]]
  out <- lm_robust(formula.m[[1]], data = s, se_type = "HC1")
  out2 <- lm_robust(formula.m[[1]], data = s, se_type = "HC1", alpha = 0.1)
  out3 <- lm_robust(formula.m[[2]], data = s, se_type = "HC1")
  out4 <- lm_robust(formula.m[[2]], data = s, se_type = "HC1", alpha = 0.1)
}
est <- array(NA, dim = c(length(treat), 5, 4))
x0 <- c(1:length(treat))+1

out <- lm_robust(formula.m[[1]], data = d.low, se_type = "HC1") # nat low edu
out2 <- lm_robust(formula.m[[1]], data = d.low, se_type = "HC1", alpha = 0.1) # nat low edu p-value = 0.1 
out3 <- lm_robust(formula.m[[1]], data = d.high, se_type = "HC1")
out4 <- lm_robust(formula.m[[1]], data = d.high, se_type = "HC1", alpha = 0.1)
out5 <- lm_robust(formula.m[[2]], data = d.low, se_type = "HC1") # pat low edu
out6 <- lm_robust(formula.m[[2]], data = d.low, se_type = "HC1", alpha = 0.1)
out7 <- lm_robust(formula.m[[2]], data = d.high, se_type = "HC1")
out8 <- lm_robust(formula.m[[2]], data = d.high, se_type = "HC1", alpha = 0.1)

est[,,1] <- cbind(out$coefficients[x0], out$conf.low[x0], out$conf.high[x0], out2$conf.low[x0], out2$conf.high[x0])
est[,,2] <- cbind(out5$coefficients[x0], out5$conf.low[x0], out5$conf.high[x0],out6$conf.low[x0], out6$conf.high[x0])
est[,,3] <- cbind(out3$coefficients[x0], out3$conf.low[x0], out3$conf.high[x0],out4$conf.low[x0], out4$conf.high[x0])
est[,,4] <- cbind(out7$coefficients[x0], out7$conf.low[x0], out7$conf.high[x0],out8$conf.low[x0], out8$conf.high[x0])

## Plotting Figure 4 left side - victim sentiment by education
par(mar = c(3,3,3,1))
plot(1, type = "n", xlim = c(0.5, 2.5), 
     ylim = c(-0.4,0.6), axes = FALSE, xlab = "", ylab = "")
box(); axis(2); 
axis(1, at = 1:2, labels = c("Below College","College or Above"), line = -2.2, tick = FALSE)
axis(1, at = 1.5, labels = c("Victim Sentiments"), line = -.2, tick = FALSE)
abline(h = seq(-1,1,0.1), col = "#AAAAAA50")
abline(v = seq(0,4,0.1), col = "#AAAAAA50")
abline(h = 0, lwd = 5, col = "#DDDDDD")
for (i in 1:2) {
  for (j in 1:length(treat)) {
    x0 <- i + adj[j]
    lines(c(x0,x0), c(est[j,2,i], est[j,3,i]), lwd = 1, col = "grey30")
    lines(c(x0,x0), c(est[j,4,i], est[j,5,i]), lwd = 3, col = "grey30")
    points(x0, est[j,1,i], col = "grey30", pch = mypch[j], cex = 1.5)
  }
}
legend("top",legend=tr.labs, border=NA,bty="n",ncol=3, col="grey30",
       lwd=2,pch=mypch,cex=0.8,seg.len=1, text.width = 0.4)
dev.copy(pdf, width = 8, height = 4, paste0("graphs/","victim_edu.pdf"))
dev.off()

# Plotting Figure 4 right side - patriotic sentiment by education
par(mar = c(3,3,3,1))
plot(1, type = "n", xlim = c(0.5, 2.5), 
     ylim = c(-0.4,0.6), axes = FALSE, xlab = "", ylab = "")
box(); axis(2); 
axis(1, at = 1:2, labels = c("Below College","College or Above"), line = -2.2, tick = FALSE)
axis(1, at = 1.5, labels = c("Patriotic Sentiments"), line = -.2, tick = FALSE)
abline(h = seq(-1,1,0.1), col = "#AAAAAA50")
abline(v = seq(0,4,0.1), col = "#AAAAAA50")
abline(h = 0, lwd = 5, col = "#DDDDDD")
for (i in 3:4) {
  for (j in 1:length(treat)) {
    x0 <- i + adj[j] -2 
    lines(c(x0,x0), c(est[j,2,i], est[j,3,i]), lwd = 1, col = "grey30")
    lines(c(x0,x0), c(est[j,4,i], est[j,5,i]), lwd = 3, col = "grey30")
    points(x0, est[j,1,i], col = "grey30", pch = mypch[j], cex = 1.5)
  }
}
legend("top",legend=tr.labs, border=NA,bty="n",ncol=3, col="grey30",
       lwd=1,pch=mypch,cex=0.8,seg.len=1, text.width = 0.4)
dev.copy(pdf, width = 8, height = 4, paste0("graphs/","pat_edu.pdf"))
dev.off()


# Figure 5 Treatment Effects on Perception of Intention and Support for Hawkish Foreign Policy
# estimates of perception of intention 
formula <- as.formula(paste0(outcomes[4],"~", paste(treat, collapse = "+"),"+",paste(covars, collapse = "+")))
out <- lm_robust(formula, data = d0, se_type = "HC1")
out2 <- lm_robust(formula, data = d0, se_type = "HC1", alpha = 0.1)
x0 <- c(1:length(treat))+1
est <- array(NA, dim = c(length(treat), 5, 2))
est[,,1] <- cbind(out$coefficients[x0], out$conf.low[x0], out$conf.high[x0],out2$conf.low[x0], out2$conf.high[x0])

# estimates of support for hawkish policies
formula2 <- as.formula(paste0(outcomes[5],"~",paste(treat, collapse = "+"),"+",paste(covars, collapse = "+")))
out3 <- lm_robust(formula2, data = d0, se_type = "HC1")
out4 <- lm_robust(formula2, data = d0, se_type = "HC1", alpha = 0.1)
x0 <- c(1:length(treat))+1
est[,,2] <- cbind(out3$coefficients[x0], out3$conf.low[x0], out3$conf.high[x0],out4$conf.low[x0], out4$conf.high[x0])

## plotting
pdf(paste0("graphs/","perc_hawk.pdf"), width = 8, height = 4)
par(mar = c(3,3,3,1))
plot(1, type = "n", xlim = c(0.5, 2.5), ylim = c(-0.4,0.6), axes = FALSE)
box(); axis(2); 
axis(1, at = 1:2, labels = y.labs[4:5], line = -.2, tick = FALSE)
abline(h = seq(-1,1,0.1), col = "#AAAAAA50") # treatment line color
abline(v = seq(0,4,0.1), col = "#AAAAAA50")
abline(h = 0, lwd = 5, col = "#DDDDDD")
for (i in 1:2) {
  for (j in 1:length(treat)) {
    x0 <- i + adj[j]
    lines(c(x0,x0), c(est[j,2,i], est[j,3,i]), lwd = 1, col = "grey30")
    lines(c(x0,x0), c(est[j,4,i], est[j,5,i]), lwd = 3, col = "grey30")
    points(x0, est[j,1,i], col = "grey30", pch = mypch[j], cex = 1.5)
  }
}
legend("top",legend=tr.labs, border=NA,bty="n",ncol=3, col = "grey30",
       lwd=2,pch=mypch,cex=1,seg.len=2)
graphics.off()


# Figure 6 Effects on Perception of Intention and Support for Hawkish Policies by Education
d.high <- d0[which(d0$eduyr>=16),]
d.low <- d0[which(d0$eduyr<16),]
d.list <- list(d.low, d.high)
sample.labs <- c("Below College","College or Above")

# estimates of perception of intention and support for hawkish policies
for (i in 1:2) {
  d0[,outcomes[i]] <- scale(d0[,outcomes[i]])
  d.list <- list(d.low, d.high)
  formula.m <- c(as.formula(paste0(outcomes[4],"~", paste(treat, collapse = "+"),"+",paste(covars, collapse = "+"))), 
                 as.formula(paste0(outcomes[5],"~", paste(treat, collapse = "+"),"+",paste(covars, collapse = "+"))))
}
est <- array(NA, dim = c(length(treat), 5, 4))
x0 <- c(1:length(treat))+1
out <- lm_robust(formula.m[[1]], data = d.low, se_type = "HC1")
out2 <- lm_robust(formula.m[[1]], data = d.low, se_type = "HC1", alpha = 0.1)
out3 <- lm_robust(formula.m[[1]], data = d.high, se_type = "HC1")
out4 <- lm_robust(formula.m[[1]], data = d.high, se_type = "HC1", alpha = 0.1)
out5 <- lm_robust(formula.m[[2]], data = d.low, se_type = "HC1")
out6 <- lm_robust(formula.m[[2]], data = d.low, se_type = "HC1", alpha = 0.1)
out7 <- lm_robust(formula.m[[2]], data = d.high, se_type = "HC1")
out8 <- lm_robust(formula.m[[2]], data = d.high, se_type = "HC1", alpha = 0.1)
est[,,1] <- cbind(out$coefficients[x0], out$conf.low[x0], out$conf.high[x0], out2$conf.low[x0], out2$conf.high[x0])
est[,,2] <- cbind(out3$coefficients[x0], out3$conf.low[x0], out3$conf.high[x0],out4$conf.low[x0], out4$conf.high[x0])
est[,,3] <- cbind(out5$coefficients[x0], out5$conf.low[x0], out5$conf.high[x0],out6$conf.low[x0], out6$conf.high[x0])
est[,,4] <- cbind(out7$coefficients[x0], out7$conf.low[x0], out7$conf.high[x0],out8$conf.low[x0], out8$conf.high[x0])
## plotting Figure 6 left side - perception of intention
par(mar = c(3,3,3,1))
plot(1, type = "n", xlim = c(0.5, 2.5), 
     ylim = c(-0.4,0.6), axes = FALSE, xlab = "", ylab = "")
box(); axis(2); 
axis(1, at = 1:2, labels = c("Below College","College or Above"), line = -2.2, tick = FALSE)
axis(1, at = 1.5, labels = c("Perception of Intention"), line = -.2, tick = FALSE)
abline(h = seq(-1,1,0.1), col = "#AAAAAA50")
abline(v = seq(0,4,0.1), col = "#AAAAAA50")
abline(h = 0, lwd = 5, col = "#DDDDDD")
for (i in 1:2) {
  for (j in 1:length(treat)) {
    x0 <- i + adj[j]
    lines(c(x0,x0), c(est[j,2,i], est[j,3,i]), lwd = 1, col = "grey30")
    lines(c(x0,x0), c(est[j,4,i], est[j,5,i]), lwd = 3, col = "grey30")
    points(x0, est[j,1,i], col = "grey30", pch = mypch[j], cex = 1.5)
  }
}
legend("top",legend=tr.labs, border=NA,bty="n",ncol=3, col="grey30",
       lwd=2,pch=mypch,cex=0.8,seg.len=1)
dev.copy(pdf, width = 8, height = 4, paste0("graphs/","perc_edu.pdf"))
dev.off()

## plotting Figure 6 right side - support for hawkish policies
par(mar = c(3,3,3,1))
plot(1, type = "n", xlim = c(0.5, 2.5), 
     ylim = c(-0.4,0.6), axes = FALSE, xlab = "", ylab = "")
box(); axis(2); 
axis(1, at = 1:2, labels = c("Below College","College or Above"), line = -2.2, tick = FALSE)
axis(1, at = 1.5, labels = c("Support for Hawkish Policies"), line = -.2, tick = FALSE)
abline(h = seq(-1,1,0.1), col = "#AAAAAA50")
abline(v = seq(0,4,0.1), col = "#AAAAAA50")
abline(h = 0, lwd = 5, col = "#DDDDDD")
for (i in 3:4) {
  for (j in 1:length(treat)) {
    x0 <- i + adj[j] -2 
    lines(c(x0,x0), c(est[j,2,i], est[j,3,i]), lwd = 1, col = "grey30")
    lines(c(x0,x0), c(est[j,4,i], est[j,5,i]), lwd = 3, col = "grey30")
    points(x0, est[j,1,i], col = "grey30", pch = mypch[j], cex = 1.5)
  }
}
legend("top",legend=tr.labs, border=NA,bty="n",ncol=3, col="grey30",
       lwd=2,pch=mypch,cex=0.8,seg.len=1)
dev.copy(pdf, width = 8, height = 4, paste0("graphs/","hawk_edu.pdf"))
dev.off()


# Figure 7 Treatment Effects on Support for the Political System
d.high <- d0[which(d0$eduyr>=16),]
d.low <- d0[which(d0$eduyr<16),]
d.list <- list(d0, d.low, d.high)
sample.labs <- c("Full Sample","Below College","College or Above")

# estimates of support for the political system
est <- array(NA, dim = c(length(treat), 5, 3))
for (m in 1:length(d.list)) {
  formula.m <- as.formula(paste0(outcomes[3],"~",paste(treat, collapse = "+"),"+",paste(covars, collapse = "+")))
  s <- d.list[[m]]
  out <- lm_robust(formula.m, data = s, se_type = "HC1")
  out2 <- lm_robust(formula.m, data = s, se_type = "HC1", alpha = 0.1)
  x0 <- c(1:length(treat))+1
  est[,,m] <- cbind(out$coefficients[x0], out$conf.low[x0], out$conf.high[x0], out2$conf.low[x0], out2$conf.high[x0])
}
## plotting
pdf(paste0("graphs/","sys_edu.pdf"), width = 8, height = 4)
par(mar = c(3,3,3,1))
plot(1, type = "n", xlim = c(0.5, 3.5), ylim = c(-0.5,0.6), axes = FALSE)
box(); axis(2); 
axis(1, at = 1:3, labels = sample.labs, line = -.2, tick = FALSE)
abline(h = seq(-1,1,0.1), col = "#AAAAAA50")
abline(v = seq(0,4,0.1), col = "#AAAAAA50")
abline(h = 0, lwd = 5, col = "#DDDDDD")
for (i in 1:3) {
  for (j in 1:length(treat)) {
    x0 <- i + adj[j]
    lines(c(x0,x0), c(est[j,2,i], est[j,3,i]), lwd = 1, col = "grey30")
    lines(c(x0,x0), c(est[j,4,i], est[j,5,i]), lwd = 3, col = "grey30")
    points(x0, est[j,1,i], col = "grey30", pch = mypch[j], cex = 1.5)
  }
}
legend("top",legend=tr.labs, border=NA,bty="n",ncol=3, col="grey30",
       lwd=2,pch=mypch,cex=1,seg.len=2)
graphics.off()

#############################
# A.8 Mediation Analysis
#############################


# victim sentiment + hawkish
nat.fit <- lm(as.formula(paste0(outcomes[1],"~", paste(treat, collapse = "+"),"+",paste(covars, collapse = "+"))), data = d0)
nat.hawk.fit <- lm(hawk ~ nat + tr_victim+tr_liberal+tr_olympic+ideology_nati+ideology_poli+ideology_econ+age+age_sq+female+hschool+jcollege+college+factor(prov)+income+substatus+knwl+minority+ccp,
                   data = d0)
nat.hawk.out <- mediate(nat.fit, nat.hawk.fit, treat = "tr_victim", mediator = "nat", robustSE = TRUE, sims = 100)
pdf(paste0("graphs/","nat_hawk.pdf"), width = 4, height = 4)
plot(nat.hawk.out, main = "Support for Hawkish Foreign Policy")
graphics.off()

# victim sentiment + perception
nat.perc.fit <- lm(perception ~ nat + tr_victim+tr_liberal+tr_olympic+ideology_nati+ideology_poli+ideology_econ+age+age_sq+female+hschool+jcollege+college+factor(prov)+income+substatus+knwl+minority+ccp,
                   data = d0)
nat.perc.out <- mediate(nat.fit, nat.perc.fit, treat = "tr_victim", mediator = "nat", robustSE = TRUE, sims = 100)
pdf(paste0("graphs/","nat_perc.pdf"),  width = 4, height = 4)
plot(nat.perc.out, main = "Perception of Intention")
graphics.off()

# pat hawkish
pat.fit <- lm(as.formula(paste0(outcomes[2],"~",paste(treat, collapse = "+"),"+",paste(covars, collapse = "+"))), data = d0)
pat.hawk.fit <- lm(hawk ~ pat + tr_victim+tr_liberal+tr_olympic+ideology_nati+ideology_poli+ideology_econ+age+age_sq+female+hschool+jcollege+college+factor(prov)+income+substatus+knwl+minority+ccp,
                   data = d0)
pat.hawk.out <- mediate(pat.fit, pat.hawk.fit, treat = "tr_victim", mediator = "pat", robustSE = TRUE, sims = 100)
pdf(paste0("graphs/","pat_hawk.pdf"), width = 4, height = 4)
plot(pat.hawk.out, main = "Support for Hawkish Foreign Policy")
graphics.off()

# pat perception
pat.perc.fit <- lm(perception ~ pat + tr_victim+tr_liberal+tr_olympic+ideology_nati+ideology_poli+ideology_econ+age+age_sq+female+hschool+jcollege+college+factor(prov)+income+substatus+knwl+minority+ccp,
                   data = d0)
pat.perc.out <- mediate(pat.fit, pat.perc.fit, treat = "tr_victim", mediator = "pat", robustSE = TRUE, sims = 100)
pdf(paste0("graphs/","pat_perc.pdf"), width = 4, height = 4)
plot(pat.perc.out,main = "Perception of Intention")
graphics.off()